rm(list=ls())
library(arfima)
library(matrixStats)
library(lattice)
#
# Create Function
#

phi <- .2
theta <- .21
d <- .2
i <- 0
x <- 100
nsim <-10

experiment <- function(phi, theta, d, i, x,nsim)
{
  res <- vector(mode="list", length=nsim)
  for(m in 1:nsim){
    # Simulate
    ts <- arfima.sim(x,model = list(phi = phi, 
                                    theta=theta, 
                                    dfrac=d, 
                                    dint =i))
    # Fit and Summarize
    arfm <- summary(arfima(ts,order = c(1, 0, 1),dmean=FALSE))
    arfm2 <- arfm[5]
    arfm3 <- arfm2$coef[[1]]
    
    #
    res[[m]] <- c(mod = arfm2,d)
  }
  class(res) <- "arfsim"
  return(res)
}


#
# Declare values
#

phi <- c(.5)
theta <- c(.3)
d <- c(0,.1,.2,.3,.4,.45)
i <- c(0)
x <- c(50,100,250,500,1000,1500)
nsim <- 100

eg <- expand.grid(x,phi,theta,d,i)
eg  
colnames(eg) <-c("x","phi","theta","d","i")

allexp <- vector(mode="list", length=nrow(eg))

#
# Run Experiments
#


for(v in 1:nrow(eg)){
  allexp[[v]] <- experiment(x =eg[v,1],
                            phi = eg[v,2], 
                            theta=eg[v,3], 
                            d=eg[v,4], 
                            i =eg[v,5],
                            nsim=nsim)
}








#
# Summarize Results
#

# Output in a matrix (eg,sim)

output<- as.matrix(allexp)
output.2 <- do.call(rbind,allexp)
output.3 <- apply(output.2,c(1,2),function(x)x[[1]]$mod.coef[1])
output.4 <- matrix(lapply(output.3,function(x)x[[1]]),nrow=nrow(eg),ncol=nsim)

output.4[[1]]

matrix.index <- expand.grid(x,d)
matrix.index

# Matrix of D estimates
output.5 <- matrix(as.numeric(lapply(output.4,function(x)x[3])),nrow=nrow(eg),ncol=nsim)
rownames(output.5) 
colnames(output.5) <- (1:nsim)

# Matrix of D Standard Errors
output.6 <- matrix(as.numeric(lapply(output.4,function(x)x[7])),nrow=nrow(eg),ncol=nsim)
rownames(output.6) 
colnames(output.6) <- (1:nsim)

# Matrix of D
output.d <- apply(output.2,c(1,2),function(x)x[[1]][[2]])
rownames(output.d) 
colnames(output.d) <- (1:nsim)



#################
## D Estimates ##
#################

#d=.0#

d.50.0 <- output.5[1,]                  
d.100.0 <- output.5[2,]  
d.250.0 <- output.5[3,]
d.500.0 <- output.5[4,]
d.1000.0 <- output.5[5,]
d.1500.0 <- output.5[6,]

e.50.0 <- output.6[1,]                  
e.100.0 <- output.6[2,]  
e.250.0 <- output.6[3,]
e.500.0 <- output.6[4,]
e.1000.0 <- output.6[5,]
e.1500.0 <- output.6[6,]

#d=.1#

d.50.1 <- output.5[7,]                  
d.100.1 <- output.5[8,]  
d.250.1 <- output.5[9,]
d.500.1 <- output.5[10,]
d.1000.1 <- output.5[11,]
d.1500.1 <- output.5[12,]

e.50.1 <- output.6[7,]                  
e.100.1 <- output.6[8,]  
e.250.1 <- output.6[9,]
e.500.1 <- output.6[10,]
e.1000.1 <- output.6[11,]
e.1500.1 <- output.6[12,]


#d=.2#

d.50.2 <- output.5[13,]                  
d.100.2 <- output.5[14,]  
d.250.2 <- output.5[15,]
d.500.2 <- output.5[16,]
d.1000.2 <- output.5[17,]
d.1500.2 <- output.5[18,]

e.50.2 <- output.6[13,]                  
e.100.2 <- output.6[14,]  
e.250.2 <- output.6[15,]
e.500.2 <- output.6[16,]
e.1000.2 <- output.6[17,]
e.1500.2 <- output.6[18,]


#d=.3#

d.50.3 <- output.5[19,]                  
d.100.3 <- output.5[20,]  
d.250.3 <- output.5[21,]
d.500.3 <- output.5[22,]
d.1000.3 <- output.5[23,]
d.1500.3 <- output.5[24,]

e.50.3 <- output.6[19,]                  
e.100.3 <- output.6[20,]  
e.250.3 <- output.6[21,]
e.500.3 <- output.6[22,]
e.1000.3 <- output.6[23,]
e.1500.3 <- output.6[24,]


#d=.4#

d.50.4 <- output.5[25,]                  
d.100.4 <- output.5[26,]  
d.250.4 <- output.5[27,]
d.500.4 <- output.5[28,]
d.1000.4 <- output.5[29,]
d.1500.4 <- output.5[30,]

e.50.4 <- output.6[25,]                  
e.100.4 <- output.6[26,]  
e.250.4 <- output.6[27,]
e.500.4 <- output.6[28,]
e.1000.4 <- output.6[29,]
e.1500.4 <- output.6[30,]

#d=.45#


d.50.45 <- output.5[31,]                  
d.100.45 <- output.5[32,]  
d.250.45 <- output.5[33,]
d.500.45 <- output.5[34,]
d.1000.45 <- output.5[35,]
d.1500.45 <- output.5[36,]

e.50.45 <- output.6[31,]                  
e.100.45 <- output.6[32,]  
e.250.45 <- output.6[33,]
e.500.45 <- output.6[34,]
e.1000.45 <- output.6[35,]
e.1500.45 <- output.6[36,]

#Average Estimates
rowMeans(output.5)

#################
## Average CIs ##
#################

#
# d.0
#


#50

d.low.50.0 <- c(d.50.0 - 2.009*(e.50.0))
d.high.50.0 <- c(d.50.0 + 2.009*(e.50.0))
d.is.50.0 <- cbind(d.low.50.0,d.50.0,d.high.50.0)
d.is.50.0
l.1.50.0 <- colMeans(d.is.50.0)

#100

d.low.100.0 <- c(d.100.0 - 1.984*(e.100.0))
d.high.100.0 <- c(d.100.0 + 1.984*(e.100.0))
d.is.100.0 <- cbind(d.low.100.0,d.100.0,d.high.100.0)
d.is.100.0
l.2.100.0 <- colMeans(d.is.100.0)

#250

d.low.250.0 <- c(d.250.0 - 1.96*(e.250.0))
d.high.250.0 <- c(d.250.0 + 1.96*(e.250.0))
d.is.250.0 <- cbind(d.low.250.0,d.250.0,d.high.250.0)
d.is.250.0
l.3.250.0 <- colMeans(d.is.250.0)

#500

d.low.500.0 <- c(d.500.0 - 1.96*(e.500.0))
d.high.500.0 <- c(d.500.0 + 1.96*(e.500.0))
d.is.500.0 <- cbind(d.low.500.0,d.500.0,d.high.500.0)
d.is.500.0
l.4.500.0 <- colMeans(d.is.500.0)


#1,000

d.low.1000.0 <- c(d.1000.0 - 1.96*(e.1000.0))
d.high.1000.0 <- c(d.1000.0 + 1.96*(e.1000.0))
d.is.1000.0 <- cbind(d.low.1000.0,d.1000.0,d.high.1000.0)
d.is.1000.0
l.5.1000.0 <- colMeans(d.is.1000.0)


#1,500

d.low.1500.0 <- c(d.1500.0 - 1.96*(e.1500.0))
d.high.1500.0 <- c(d.1500.0 + 1.96*(e.1500.0))
d.is.1500.0 <- cbind(d.low.1500.0,d.1500.0,d.high.1500.0)
d.is.1500.0
l.6.1500.0 <- colMeans(d.is.1500.0)


t.o.0 <- cbind(l.1.50.0,l.2.100.0,l.3.250.0,l.4.500.0,l.5.1000.0,l.6.1500.0)
table.out.0 <- t(t.o.0)
rownames(table.out.0) 
colnames(table.out.0) <- c("d.low","d.est","d.high")
table.out.0

# Plots of Distributions
par(mfrow=c(1,1))
plot (density(d.50.0), col="red",ylim=c(0,7),xlim=c(-1,1)) 
lines (density(d.100.0), col="green") 
lines (density(d.250.0), col="blue")
lines (density(d.500.0), col="orange")
lines (density(d.1000.0), col="black")
lines (density(d.1500.0), col="grey")
abline(v=0)


#
# d.1
#


#50

d.low.50.1 <- c(d.50.1 - 2.009*(e.50.1))
d.high.50.1 <- c(d.50.1 + 2.009*(e.50.1))
d.is.50.1 <- cbind(d.low.50.1,d.50.1,d.high.50.1)
d.is.50.1
l.1.50.1 <- colMeans(d.is.50.1)

#100

d.low.100.1 <- c(d.100.1 - 1.984*(e.100.1))
d.high.100.1 <- c(d.100.1 + 1.984*(e.100.1))
d.is.100.1 <- cbind(d.low.100.1,d.100.1,d.high.100.1)
d.is.100.1
l.2.100.1 <- colMeans(d.is.100.1)

#250

d.low.250.1 <- c(d.250.1 - 1.96*(e.250.1))
d.high.250.1 <- c(d.250.1 + 1.96*(e.250.1))
d.is.250.1 <- cbind(d.low.250.1,d.250.1,d.high.250.1)
d.is.250.1
l.3.250.1 <- colMeans(d.is.250.1)

#500

d.low.500.1 <- c(d.500.1 - 1.96*(e.500.1))
d.high.500.1 <- c(d.500.1 + 1.96*(e.500.1))
d.is.500.1 <- cbind(d.low.500.1,d.500.1,d.high.500.1)
d.is.500.1
l.4.500.1 <- colMeans(d.is.500.1)


#1,000

d.low.1000.1 <- c(d.1000.1 - 1.96*(e.1000.1))
d.high.1000.1 <- c(d.1000.1 + 1.96*(e.1000.1))
d.is.1000.1 <- cbind(d.low.1000.1,d.1000.1,d.high.1000.1)
d.is.1000.1
l.5.1000.1 <- colMeans(d.is.1000.1)


#1,500

d.low.1500.1 <- c(d.1500.1 - 1.96*(e.1500.1))
d.high.1500.1 <- c(d.1500.1 + 1.96*(e.1500.1))
d.is.1500.1 <- cbind(d.low.1500.1,d.1500.1,d.high.1500.1)
d.is.1500.1
l.6.1500.1 <- colMeans(d.is.1500.1)


t.o.1 <- cbind(l.1.50.1,l.2.100.1,l.3.250.1,l.4.500.1,l.5.1000.1,l.6.1500.1)
table.out.1 <- t(t.o.1)
rownames(table.out.1) 
colnames(table.out.1) <- c("d.low","d.est","d.high")
table.out.1

# Plots of Distributions

plot (density(d.50.1), col="red",ylim=c(0,6),xlim=c(-1,1)) 
lines (density(d.100.1), col="green") 
lines (density(d.250.1), col="blue")
lines (density(d.500.1), col="orange")
lines (density(d.1000.1), col="black")
lines (density(d.1500.1), col="grey")
abline(v=.1)


#
# d.2
#


#50

d.low.50.2 <- c(d.50.2 - 2.009*(e.50.2))
d.high.50.2 <- c(d.50.2 + 2.009*(e.50.2))
d.is.50.2 <- cbind(d.low.50.2,d.50.2,d.high.50.2)
d.is.50.2
l.1.50.2 <- colMeans(d.is.50.2)

#100

d.low.100.2 <- c(d.100.2 - 1.984*(e.100.2))
d.high.100.2 <- c(d.100.2 + 1.984*(e.100.2))
d.is.100.2 <- cbind(d.low.100.2,d.100.2,d.high.100.2)
d.is.100.2
l.2.100.2 <- colMeans(d.is.100.2)

#250

d.low.250.2 <- c(d.250.2 - 1.96*(e.250.2))
d.high.250.2 <- c(d.250.2 + 1.96*(e.250.2))
d.is.250.2 <- cbind(d.low.250.2,d.250.2,d.high.250.2)
d.is.250.2
l.3.250.2 <- colMeans(d.is.250.2)

#500

d.low.500.2 <- c(d.500.2 - 1.96*(e.500.2))
d.high.500.2 <- c(d.500.2 + 1.96*(e.500.2))
d.is.500.2 <- cbind(d.low.500.2,d.500.2,d.high.500.2)
d.is.500.2
l.4.500.2 <- colMeans(d.is.500.2)


#1,000

d.low.1000.2 <- c(d.1000.2 - 1.96*(e.1000.2))
d.high.1000.2 <- c(d.1000.2 + 1.96*(e.1000.2))
d.is.1000.2 <- cbind(d.low.1000.2,d.1000.2,d.high.1000.2)
d.is.1000.2
l.5.1000.2 <- colMeans(d.is.1000.2)


#1,500

d.low.1500.2 <- c(d.1500.2 - 1.96*(e.1500.2))
d.high.1500.2 <- c(d.1500.2 + 1.96*(e.1500.2))
d.is.1500.2 <- cbind(d.low.1500.2,d.1500.2,d.high.1500.2)
d.is.1500.2
l.6.1500.2 <- colMeans(d.is.1500.2)


t.o.2 <- cbind(l.1.50.2,l.2.100.2,l.3.250.2,l.4.500.2,l.5.1000.2,l.6.1500.2)
table.out.2 <- t(t.o.2)
rownames(table.out.2) 
colnames(table.out.2) <- c("d.low","d.est","d.high")
table.out.2

# Plots of Distributions

plot (density(d.50.2), col="red",ylim=c(0,6),xlim=c(-1,1)) 
lines (density(d.100.2), col="green") 
lines (density(d.250.2), col="blue")
lines (density(d.500.2), col="orange")
lines (density(d.1000.2), col="black")
lines (density(d.1500.2), col="grey")
abline(v=.2)


#
# d.3
#


#50

d.low.50.3 <- c(d.50.3 - 2.009*(e.50.3))
d.high.50.3 <- c(d.50.3 + 2.009*(e.50.3))
d.is.50.3 <- cbind(d.low.50.3,d.50.3,d.high.50.3)
d.is.50.3
l.1.50.3 <- colMeans(d.is.50.3)

#100

d.low.100.3 <- c(d.100.3 - 1.984*(e.100.3))
d.high.100.3 <- c(d.100.3 + 1.984*(e.100.3))
d.is.100.3 <- cbind(d.low.100.3,d.100.3,d.high.100.3)
d.is.100.3
l.2.100.3 <- colMeans(d.is.100.3)

#250

d.low.250.3 <- c(d.250.3 - 1.96*(e.250.3))
d.high.250.3 <- c(d.250.3 + 1.96*(e.250.3))
d.is.250.3 <- cbind(d.low.250.3,d.250.3,d.high.250.3)
d.is.250.3
l.3.250.3 <- colMeans(d.is.250.3)

#500

d.low.500.3 <- c(d.500.3 - 1.96*(e.500.3))
d.high.500.3 <- c(d.500.3 + 1.96*(e.500.3))
d.is.500.3 <- cbind(d.low.500.3,d.500.3,d.high.500.3)
d.is.500.3
l.4.500.3 <- colMeans(d.is.500.3)


#1,000

d.low.1000.3 <- c(d.1000.3 - 1.96*(e.1000.3))
d.high.1000.3 <- c(d.1000.3 + 1.96*(e.1000.3))
d.is.1000.3 <- cbind(d.low.1000.3,d.1000.3,d.high.1000.3)
d.is.1000.3
l.5.1000.3 <- colMeans(d.is.1000.3)


#1,500

d.low.1500.3 <- c(d.1500.3 - 1.96*(e.1500.3))
d.high.1500.3 <- c(d.1500.3 + 1.96*(e.1500.3))
d.is.1500.3 <- cbind(d.low.1500.3,d.1500.3,d.high.1500.3)
d.is.1500.3
l.6.1500.3 <- colMeans(d.is.1500.3)


t.o.3 <- cbind(l.1.50.3,l.2.100.3,l.3.250.3,l.4.500.3,l.5.1000.3,l.6.1500.3)
table.out.3 <- t(t.o.3)
rownames(table.out.3) 
colnames(table.out.3) <- c("d.low","d.est","d.high")
table.out.3

# Plots of Distributions

plot (density(d.50.3), col="red",ylim=c(0,2.5),xlim=c(-1,1)) 
lines (density(d.100.3), col="green") 
lines (density(d.250.3), col="blue")
lines (density(d.500.3), col="orange")
lines (density(d.1000.3), col="black")
lines (density(d.1500.3), col="grey")
abline(v=.3)


#
# d.4
#


#50

d.low.50.4 <- c(d.50.4 - 2.009*(e.50.4))
d.high.50.4 <- c(d.50.4 + 2.009*(e.50.4))
d.is.50.4 <- cbind(d.low.50.4,d.50.4,d.high.50.4)
d.is.50.4
l.1.50.4 <- colMeans(d.is.50.4)

#100

d.low.100.4 <- c(d.100.4 - 1.984*(e.100.4))
d.high.100.4 <- c(d.100.4 + 1.984*(e.100.4))
d.is.100.4 <- cbind(d.low.100.4,d.100.4,d.high.100.4)
d.is.100.4
l.2.100.4 <- colMeans(d.is.100.4)

#250

d.low.250.4 <- c(d.250.4 - 1.96*(e.250.4))
d.high.250.4 <- c(d.250.4 + 1.96*(e.250.4))
d.is.250.4 <- cbind(d.low.250.4,d.250.4,d.high.250.4)
d.is.250.4
l.3.250.4 <- colMeans(d.is.250.4)

#500

d.low.500.4 <- c(d.500.4 - 1.96*(e.500.4))
d.high.500.4 <- c(d.500.4 + 1.96*(e.500.4))
d.is.500.4 <- cbind(d.low.500.4,d.500.4,d.high.500.4)
d.is.500.4
l.4.500.4 <- colMeans(d.is.500.4)


#1,000

d.low.1000.4 <- c(d.1000.4 - 1.96*(e.1000.4))
d.high.1000.4 <- c(d.1000.4 + 1.96*(e.1000.4))
d.is.1000.4 <- cbind(d.low.1000.4,d.1000.4,d.high.1000.4)
d.is.1000.4
l.5.1000.4 <- colMeans(d.is.1000.4)


#1,500

d.low.1500.4 <- c(d.1500.4 - 1.96*(e.1500.4))
d.high.1500.4 <- c(d.1500.4 + 1.96*(e.1500.4))
d.is.1500.4 <- cbind(d.low.1500.4,d.1500.4,d.high.1500.4)
d.is.1500.4
l.6.1500.4 <- colMeans(d.is.1500.4)


t.o.4 <- cbind(l.1.50.4,l.2.100.4,l.3.250.4,l.4.500.4,l.5.1000.4,l.6.1500.4)
table.out.4 <- t(t.o.4)
rownames(table.out.4) 
colnames(table.out.4) <- c("d.low","d.est","d.high")
table.out.4

# Plots of Distributions

plot (density(d.50.4), col="red",ylim=c(0,3),xlim=c(-1,1)) 
lines (density(d.100.4), col="green") 
lines (density(d.250.4), col="blue")
lines (density(d.500.4), col="orange")
lines (density(d.1000.4), col="black")
lines (density(d.1500.4), col="grey")
abline(v=.4)



#
# d.45
#


#50

d.low.50.45 <- c(d.50.45 - 2.009*(e.50.45))
d.high.50.45 <- c(d.50.45 + 2.009*(e.50.45))
d.is.50.45 <- cbind(d.low.50.45,d.50.45,d.high.50.45)
d.is.50.45
l.1.50.45 <- colMeans(d.is.50.45)

#100

d.low.100.45 <- c(d.100.45 - 1.984*(e.100.45))
d.high.100.45 <- c(d.100.45 + 1.984*(e.100.45))
d.is.100.45 <- cbind(d.low.100.45,d.100.45,d.high.100.45)
d.is.100.45
l.2.100.45 <- colMeans(d.is.100.45)

#250

d.low.250.45 <- c(d.250.45 - 1.96*(e.250.45))
d.high.250.45 <- c(d.250.45 + 1.96*(e.250.45))
d.is.250.45 <- cbind(d.low.250.45,d.250.45,d.high.250.45)
d.is.250.45
l.3.250.45 <- colMeans(d.is.250.45)

#500

d.low.500.45 <- c(d.500.45 - 1.96*(e.500.45))
d.high.500.45 <- c(d.500.45 + 1.96*(e.500.45))
d.is.500.45 <- cbind(d.low.500.45,d.500.45,d.high.500.45)
d.is.500.45
l.4.500.45 <- colMeans(d.is.500.45)


#1,000

d.low.1000.45 <- c(d.1000.45 - 1.96*(e.1000.45))
d.high.1000.45 <- c(d.1000.45 + 1.96*(e.1000.45))
d.is.1000.45 <- cbind(d.low.1000.45,d.1000.45,d.high.1000.45)
d.is.1000.45
l.5.1000.45 <- colMeans(d.is.1000.45)


#1,500

d.low.1500.45 <- c(d.1500.45 - 1.96*(e.1500.45))
d.high.1500.45 <- c(d.1500.45 + 1.96*(e.1500.45))
d.is.1500.45 <- cbind(d.low.1500.45,d.1500.45,d.high.1500.45)
d.is.1500.45
l.6.1500.45 <- colMeans(d.is.1500.45)


t.o.45 <- cbind(l.1.50.45,l.2.100.45,l.3.250.45,l.4.500.45,l.5.1000.45,l.6.1500.45)
table.out.45 <- t(t.o.45)
rownames(table.out.45) 
colnames(table.out.45) <- c("d.low","d.est","d.high")
table.out.45

# Plots of Distributions

plot (density(d.50.45), col="red",ylim=c(0,3.5),xlim=c(-1,1)) 
lines (density(d.100.45), col="green") 
lines (density(d.250.45), col="blue")
lines (density(d.500.45), col="orange")
lines (density(d.1000.45), col="black")
lines (density(d.1500.45), col="grey")
abline(v=.45)



#
# Final Output ARIMA(1,d,1) phi=.5 theta=.3
#

table.out.0
table.out.1
table.out.2
table.out.3
table.out.4
table.out.45

table.out.0
summary(d.50.0)
summary(d.100.0)
summary(d.250.0)
summary(d.500.0)
summary(d.1000.0)
summary(d.1500.0)

table.out.1
summary(d.50.1)
summary(d.100.1)
summary(d.250.1)
summary(d.500.1)
summary(d.1000.1)
summary(d.1500.1)

table.out.2
summary(d.50.2)
summary(d.100.2)
summary(d.250.2)
summary(d.500.2)
summary(d.1000.2)
summary(d.1500.2)

table.out.3
summary(d.50.3)
summary(d.100.3)
summary(d.250.3)
summary(d.500.3)
summary(d.1000.3)
summary(d.1500.3)

table.out.4
summary(d.50.4)
summary(d.100.4)
summary(d.250.4)
summary(d.500.4)
summary(d.1000.4)
summary(d.1500.4)

table.out.45
summary(d.50.45)
summary(d.100.45)
summary(d.250.45)
summary(d.500.45)
summary(d.1000.45)
summary(d.1500.45)


par(mfrow=c(2,3))
plot (density(d.50.0), col="red",ylim=c(0,7),xlim=c(-1,1),main="d=0") 
lines (density(d.100.0), col="green") 
lines (density(d.250.0), col="blue")
lines (density(d.500.0), col="orange")
lines (density(d.1000.0), col="black")
lines (density(d.1500.0), col="grey")
abline(v=0)

plot (density(d.50.1), col="red",ylim=c(0,6),xlim=c(-1,1),main="d=.1") 
lines (density(d.100.1), col="green") 
lines (density(d.250.1), col="blue")
lines (density(d.500.1), col="orange")
lines (density(d.1000.1), col="black")
lines (density(d.1500.1), col="grey")
abline(v=.1)

plot (density(d.50.2), col="red",ylim=c(0,6),xlim=c(-1,1),main="d=.2") 
lines (density(d.100.2), col="green") 
lines (density(d.250.2), col="blue")
lines (density(d.500.2), col="orange")
lines (density(d.1000.2), col="black")
lines (density(d.1500.2), col="grey")
abline(v=.2)

plot (density(d.50.3), col="red",ylim=c(0,2.5),xlim=c(-1,1),main="d=.3") 
lines (density(d.100.3), col="green") 
lines (density(d.250.3), col="blue")
lines (density(d.500.3), col="orange")
lines (density(d.1000.3), col="black")
lines (density(d.1500.3), col="grey")
abline(v=.3)

plot (density(d.50.4), col="red",ylim=c(0,3),xlim=c(-1,1),main="d=.4") 
lines (density(d.100.4), col="green") 
lines (density(d.250.4), col="blue")
lines (density(d.500.4), col="orange")
lines (density(d.1000.4), col="black")
lines (density(d.1500.4), col="grey")
abline(v=.4)

plot (density(d.50.45), col="red",ylim=c(0,3.5),xlim=c(-1,1),main="d=.45") 
lines (density(d.100.45), col="green") 
lines (density(d.250.45), col="blue")
lines (density(d.500.45), col="orange")
lines (density(d.1000.45), col="black")
lines (density(d.1500.45), col="grey")
abline(v=.45)
